#include "matrix.h"
#include <cstdio>
#include <iostream>
#include <cmath>
#include <algorithm>
#pragma warning(disable : 4996)

Matrix ::Matrix(int h, int w)
{
    height = h;
    width = w;
    data = nullptr;
    da=nullptr;
}
Matrix ::Matrix()
{
    height=0;
    width=0;
    data=nullptr;
    da=nullptr;
}
Matrix::~Matrix()
{
    delete[]data;
    if(da!=nullptr)
    {
        for (int i = 0; i < height; i++)
            {
                delete[]da[i];
            }
            delete[]da;
    }
}
void Matrix::Readimage(QString  filename)
{
 /*   std::string str = filename.toStdString();
    const char* ch = str.c_str();
    FILE* fp = NULL;
    fp = fopen(ch, "rb");
    fread(&bm, sizeof(bm), 1, fp);
    fread(&bi, sizeof(bi), 1, fp);
     da = new double* [height];
      for (int i=0;i<height;i++)
     {
         da[i] = new double [width*3];
          for (int j=0;j<width;j++)
         {
              da[i][j*3]=*(data+2);
              da[i][j*3+1]=*(data+1);
              da[i][j*3+2]=*(data)  ;
              data+=4;
         }
     }*/
    std::string str = filename.toStdString();
    const char* ch = str.c_str();
    FILE* fp = NULL;
    unsigned char x=0;
    fp = fopen(ch, "rb");
    fread(&bm, sizeof(bm), 1, fp);
    fread(&bi, sizeof(bi), 1, fp);
    height=bi.bihight;
    width=bi.biwidth;
     da = new double* [height];
      for (int i=0;i<height;i++)
     {
            da[i] = new double[width * 3];
          for (int j = 0; j < (width * 3 + 3) / 4 * 4 - 2 * width; j++)

         {


             if (j < width)
             {
                 fread(&x, sizeof(unsigned char), 1, fp);
                 da[i][j * 3 + 0] = x;
                 fread(&x, sizeof(unsigned char), 1, fp);
                 da[i][j * 3 + 1] = x;
                 fread(&x, sizeof(unsigned char), 1, fp);
                 da[i][j * 3 + 2] = x;

             }
             else
             {
                 fread(&x, 1, 1, fp);
             }
         }
     }
      fclose(fp);
}
void Matrix::doubletosingle()
{
    for (int i=height-1;i>=0;i--)
       {
            for (int j=0;j<width;j++)
           {
                *(data+2)=da[i][j*3+2];
                *(data+1)=da[i][j*3+1];
                *(data )=da[i][j*3 ];
                data+=4;
           }
       }
}
void Matrix::Writeimage(QString filename)
{
       /* std::string str = filename.toStdString();
        const char* ch = str.c_str();
        FILE* fp = NULL;
        fp = fopen(ch, "wb");
        fwrite(&bm, sizeof(bm), 1, fp);
        fwrite(&bi, sizeof(bi), 1, fp);
        unsigned char x = 0;
        data+=(4*height*width-4);
        data-=(4*(width-1));
        for (int i = 0; i < height; i++)
        {

            for (int j = 0; j < (width * 3 + 3) / 4 * 4 - 2 * width; j++)
            {
                if (j < width)
                {
                    x = *(data);
                    fwrite(&x, 1, 1, fp);
                    x = *(data+1);
                    fwrite(&x, 1, 1, fp);
                    x = *(data +2);
                    fwrite(&x, 1, 1, fp);
                    data+=4;
                }
                else
                {
                    x = 0;
                    fwrite(&x, 1, 1, fp);
                }
            }
            data-=width*4;
            data-=width*4;
        }
        fclose(fp);*/
    std::string str = filename.toStdString();
    const char* ch = str.c_str();
    FILE* fp = NULL;
    fp = fopen(ch, "wb");
    fwrite(&bm, sizeof(bm), 1, fp);
    fwrite(&bi, sizeof(bi), 1, fp);
    unsigned char x = 0;

        for (int i = 0; i < height; i++)
        {

            for (int j = 0; j < (width * 3 + 3) / 4 * 4 - 2 * width; j++)
            {
                if (j < width)
                {
                    x = da[i][j * 3 + 0];
                    fwrite(&x, 1, 1, fp);
                    x = da[i][j * 3 + 1];
                    fwrite(&x, 1, 1, fp);
                    x = da[i][j * 3 + 2];
                    fwrite(&x, 1, 1, fp);
                }
                else
                {
                    x = 0;
                    fwrite(&x, 1, 1, fp);
                }
            }
        }
    fclose(fp);
}
void Matrix::copy(const Matrix &m)
{
    for(int i=0;i<height;i++)
    {
        delete []da[i];

    }
    delete da;
    height=m.height;
    width=m.width;
    bi=m.bi;
    bm=m.bm;
    da=new double *[height];
    for(int i=0;i<height;i++)
    {
        da[i]=new double [width*3];
        for(int j=0;j<width;j++)
        {
           da[i][j*3]=m.da[i][j*3];
           da[i][j*3+1]=m.da[i][j*3+1];
           da[i][j*3+2]=m.da[i][j*3+2];
        }
    }


}
void Matrix::light(double t)
{
    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j < width; j++)

        {
            da[i][j * 3 + 2] *=t;
            da[i][j * 3 + 1] *=t;
            da[i][j * 3 + 0] *=t;
            if(da[i][j * 3 + 0]>=255)
                da[i][j * 3 + 0]=255;
            if(da[i][j * 3 + 1]>=255)
                da[i][j * 3 + 1]=255;
            if(da[i][j * 3 + 2]>=255)
               da[i][j * 3 + 2]=255;
        }
    }
}
void Matrix::Flip(int code)
{
    if (code)
    {
        double* m = nullptr;
    //	for (int i = 0; i < width * 3; i++)
            for (int j = 0; j < (height + 1) / 2; j++)
            {
                m = da[height - 1 - j];
                da[height - 1 - j] = da[j];
                da[j] = m;
            }
    }
    else
    {
        double t = 0;
        for (int i = 0; i < height; i++)
            for (int j = 0; j < (width-1) / 2; j++)
            {
                int x = j * 3 + 0, y = j * 3 + 1, z = j * 3 + 2;
                t = da[i][(width - 1 - j)*3];
                da[i][(width - 1 - j)*3] = da[i][x];
                da[i][x] = t;

                t = da[i][(width  - 1 - j)*3+1];
                da[i][(width  - 1 - j)*3+1] = da[i][y];
                da[i][y] = t;

                t = da[i][(width  - 1 - j)*3+2];
                da[i][(width  - 1 - j)*3+2] = da[i][z];
                da[i][z] = t;
            }
    }
}
void Matrix::grey()
{
    double x = 0;
    for (int i = 0; i < height; i++)
    {
        for (int j = 0; j < width; j++)
        {
            x = 0;
            x = (pow(da[i][j * 3], 2.2) + pow(1.5 * da[i][j * 3 + 1], 2.2) + pow(0.6 * da[i][j * 3 + 2], 2.2)) / (1 + pow(1.5, 2.2) + pow(0.6, 2.2));
            x = pow(x, 1 / 2.2);
            da[i][j * 3]=x;
            da[i][j * 3+1]=x;
            da[i][j * 3+2]=x;
        }
    }
}
void Matrix::mosaic(int x1, int y1, int x2, int y2, int len)
{
    if(x1>x2)
    {
        int t=x1;
        x1=x2;
        x2=t;
    }
    if(y1>y2)
    {
        int t=y1;
        y1=y2;
        y2=t;
    }
        int x = x1, y = y1;
        for (; x <= x2;)
        {
            for (; y <= y2;)
            {
                if (x + len - 1 > x2)x = x2 - len + 1;
                if (y + len - 1 > y2)y = y2 - len + 1;
                for (int i = x; i < x + len; i++)
                {
                    for (int j = y; j < y + len; j++)
                    {
                        da[i][j * 3] = da[x + len - 1][(y + len - 1) * 3];
                        da[i][j * 3 + 1] = da[x + len - 1][(y + len - 1) * 3 + 1];
                        da[i][j * 3 + 2] = da[x + len - 1][(y + len - 1) * 3 + 2];
                    }
                }
                y = y + len;
            }
            x = x + len;
            y = y1;
        }

}
void Matrix::Resize(int h, int w)
{
    int i, j;
    double** data1;
    data1 = new double* [h];
    for (int i = 0; i < h; i++)
    {
        data1[i] = new double[w * 3];
    }

    int xx, yy;
    double rh, rw;
    for (i = 0; i < h; i++)
    {

        for (j = 0; j < w; j++)
        {
            xx = i * height / h; yy = j * width / w;
            rh = (double)i * height / h - (int)i * height / h;
            rw = (double)j * width / w - (int)j * width / w;
            if (xx + 1 >= height)xx = xx - 1;
            if (yy + 1 >= width)yy = yy - 1;
            if (rh == 0 && rw == 0) {
                data1[i][j * 3] = da[xx][yy * 3];
                data1[i][j * 3 + 1] = da[xx][yy * 3 + 1];
                data1[i][j * 3 + 2] = da[xx][yy * 3 + 2];
                continue;
            }
            if (rh != 0 && rw != 0) {
                data1[i][j * 3] = rh * rw * da[xx + 1][(yy + 1) * 3] + (1 - rh) * rw * da[xx][(yy + 1) * 3] + (1 - rw) * rh * da[xx + 1][yy * 3] + (1 - rh) * (1 - rw) * da[xx][yy * 3];
                data1[i][j * 3 + 1] = rh * rw * da[xx + 1][(yy + 1) * 3 + 1] + (1 - rh) * rw * da[xx][(yy + 1) * 3 + 1] + (1 - rw) * rh * da[xx + 1][yy * 3 + 1] + (1 - rh) * (1 - rw) * da[xx][yy * 3 + 1];
                data1[i][j * 3 + 2] = rh * rw * da[xx + 1][(yy + 1) * 3 + 2] + (1 - rh) * rw * da[xx][(yy + 1) * 3 + 2] + (1 - rw) * rh * da[xx + 1][yy * 3 + 1] + (1 - rh) * (1 - rw) * da[xx][yy * 3 + 2];

                continue;
            }
            if (rh != 0 && rw == 0) {
                data1[i][j * 3] = rh * da[xx + 1][yy * 3] + (1 - rh) * da[xx][yy * 3];
                data1[i][j * 3 + 1] = rh * da[xx + 1][yy * 3 + 1] + (1 - rh) * da[xx][yy * 3 + 1];
                data1[i][j * 3 + 2] = rh * da[xx + 1][yy * 3 + 2] + (1 - rh) * da[xx][yy * 3 + 2];
                continue;
            }
            if (rh == 0 && rw != 0) {
                data1[i][j * 3] = da[xx][(yy + 1) * 3] * rw + (1 - rw) * da[xx][yy * 3];
                data1[i][j * 3 + 1] = da[xx][(yy + 1) * 3 + 1] * rw + (1 - rw) * da[xx][yy * 3 + 1];
                data1[i][j * 3 + 2] = da[xx][(yy + 1) * 3 + 2] * rw + (1 - rw) * da[xx][yy * 3 + 2];
                continue;
            }

        }

    }
    for (i = 0; i < height; i++)
    {
        delete[]da[i];
    }

    delete[]da;
    da = data1;

    height = h; width = w;
    bi.bihight = height;
    bi.biwidth = width;
    bi.bisizeimage = (((width * bi.bibitcount) + 31) / 32 * 4) * height;
}
void Matrix::Cut(int x1, int y1, int x2, int y2)
{
    if(x1>x2)
    {
        int t=x1;
        x1=x2;
        x2=t;
    }
    if(y1>y2)
    {
        int t=y1;
        y1=y2;
        y2=t;
    }
        int w = x2 - x1 + 1, h = y2 - y1 + 1;
        double** data1;
        data1 = new double* [h];
        for (int i = 0; i < h; i++)
        {
            data1[i] = new double[w * 3];
            for (int j = 0; j < w; j++)
            {

                data1[i][j * 3] = da[i + x1][(j + y1) * 3];
                data1[i][j * 3 + 1] = da[i + x1][(j + y1) * 3 + 1];
                data1[i][j * 3 + 2] = da[i + x1][(j + y1) * 3 + 2];
            }
        }
        for (int i = 0; i < height; i++)
        {
            delete[]da[i];
        }
        delete[]da;

        da = data1;
        height = h; width = w;
        bi.biwidth = width; bi.bihight = height;
        bi.bisizeimage = (((width * bi.bibitcount) + 31) / 32 * 4) * height;
}
void Matrix::duibidu(double t)
{
    double x = 0, y = 0, z = 0, m = 127;
        //double r = 0, g = 0, b = 0, rb = 0, rg = 0, rr = 0;
        for (int i = 0; i < height; i++)
        {
            for (int j = 0; j < width; j++)
            {
                x = da[i][j * 3];
                y = da[i][j * 3 + 1];
                z = da[i][j * 3 + 2];
                x = x + (x - m) * t / 255;
                y = y + (y - m) * t / 255;
                z = z + (z - m) * t / 255;
                da[i][j * 3] = x;
                da[i][j * 3 + 1] = y;
                da[i][j * 3 + 2] = z;
                if (x > 256)x = 255;
                if (y > 256)y = 255;
                if (z > 256)z = 255;
                da[i][j * 3] = x;
                da[i][j * 3 + 1] = y;
                da[i][j * 3 + 2] = z;
            }
        }

}
void Matrix::Rotate(int code)
{
    int t;
        double** data1;
        //double x = 0;

        t = width;
        width = height;
        height = t;
        data1 = new double* [height];

        if (code)
        {
            for (int i = 0; i < height; i++)
            {
                data1[i] = new double[width * 3];
                for (int j = 0; j < width; j++)
                {
                    data1[i][j * 3] = da[width - 1 - j][i * 3];
                    data1[i][j * 3 + 1] = da[width - 1 - j][i * 3 + 1];
                    data1[i][j * 3 + 2] = da[width - 1 - j][i * 3 + 2];
                }

            }
        }
        else
        {
            for (int i = 0; i < height; i++)
            {
                data1[i] = new double[width * 3];
                for (int j = 0; j < width; j++)
                {
                    data1[i][j * 3] = da[j][(height - 1 - i) * 3];
                    data1[i][j * 3 + 1] = da[j][(height - 1 - i) * 3 + 1];
                    data1[i][j * 3 + 2] = da[j][(height - 1 - i) * 3 + 2];
                    if (j == width - 1)
                    {
                        //int m = 0;
                    }
                }

            }
        }
        for (int i = 0; i < width; i++)
        {
            delete[]da[i];
        }
        delete[]da;
        da = data1;


        bi.bihight = height;
        bi.biwidth = width;
        bi.bisizeimage = (((width * bi.bibitcount) + 31) / 32 * 4) * height;
}
void Matrix::gaosilvboqi(int ksize,double strength)
{
    double pi = 3.1415926;
        double window[11][11];
        double sigma = strength;
        int center = ksize / 2;
        double x2, y2;
        for (int i = 0; i < ksize; i++)
        {
            x2 = pow(i - center, 2);
            for (int j = 0; j < ksize; j++)
            {
                y2 = pow(j - center, 2);
                double g = exp(-(x2 + y2) / (2 * sigma * sigma));
                g /= 2 * pi * sigma;
                window[i][j] = g;
            }
        }
        double k = 1 / window[0][0];
        int s = 0;
        for (int i = 0; i < ksize; i++)
        {
            for (int j = 0; j < ksize; j++)
            {
                window[i][j] *= k;
                s = s + (int)window[i][j];
            }
        }
        double sum = 0;
        double sum1 = 0, sum2 = 0;
        for (int i = 0; i < height; i++)
        {
            for (int j = 0; j < width; j++)
            {
                sum = 0;
                sum1 = 0;
                sum2 = 0;
                for (int k = i - center; k <= i + center; k++)
                {
                    for (int l = j - center; l <= j + center; l++)
                    {
                        if (l >= 0 && l < width && k >= 0 && k < height)
                        {
                            sum += (da[k][l * 3] * (int)window[k - i + center][l - j + center]);
                            sum1 += (da[k][l * 3 + 1] * (int)window[k - i + center][l - j + center]);
                            sum2 += (da[k][l * 3 + 2] * (int)window[k - i + center][l - j + center]);
                        }
                    }
                }
                sum = sum / s;
                sum1 = sum1 / s;
                sum2 = sum2 / s;
                da[i][j * 3] = sum;
                da[i][j * 3 + 1] = sum1;
                da[i][j * 3 + 2] = sum2;
            }
        }

}
void Matrix::tidu()
{
    double x = 0, y = 0;
        for (int i = 0; i < height - 1; i++)
        {
            for (int j = 0; j < width - 1; j++)
            {
                x = (da[i][j * 3 + 3] - da[i + 1][3 * j] - da[i][3 * j] + da[i + 1][3 * j + 3]) / 2;
                y = (da[i][j * 3] - da[i + 1][3 * j] + da[i][3 * j + 3] - da[i + 1][3 * j + 3]) / 2;
                da[i][3 * j] = sqrt(x * x + y * y);
                da[i][3 * j + 1] = sqrt(x * x + y * y);
                da[i][3 * j + 2] = sqrt(x * x + y * y);

            }
        }
}
void Matrix::baohedu(double percent)
{
    double x = 0, y = 0, z = 0, rbgmax = 0, rbgmin = 0,
            delta = 0, value = 0, L = 0, S = 0, alpha = 0;
        for (int i = 0; i < height; i++)
        {
            for (int j = 0; j < width; j++)
            {
                x = da[i][j * 3];
                y = da[i][j * 3 + 1];
                z = da[i][j * 3 + 2];
                rbgmax = fmax(fmax(x, y), z);
                rbgmin = fmin(fmin(x, y), z);
                delta = (rbgmax - rbgmin) / 255;
                if (rbgmax != rbgmin)
                {
                    value = (rbgmax + rbgmin) / 255;
                    L = value / 2;
                    if (L < 0.5)
                    {
                        S = delta / value;
                    }
                    else
                    {
                        S = delta / (2 - value);
                    }
                    if (percent >= 0)
                    {
                        if (percent + S >= 1)
                            alpha = S;
                        else
                            alpha = 1 - percent;
                        alpha = 1 / alpha - 1;
                        x = x + (x - L * 255) * alpha;
                        y = y + (y - L * 255) * alpha;
                        z = z + (z - L * 255) * alpha;
                    }
                    else
                    {
                        alpha = percent;
                        x = L * 255 + (x - L * 255) * (1 + alpha);
                        y = L * 255 + (x - L * 255) * (1 + alpha);
                        z = L * 255 + (x - L * 255) * (1 + alpha);
                    }
                    da[i][j * 3] = x;
                    da[i][j * 3 + 1] = y;
                    da[i][j * 3 + 2] = z;
                }
            }
        }
}
void Matrix::gamma(double t)
{
    double x = 0, y = 0, z = 0;
        t = 1 / t;
        for (int i = 0; i < height; i++)
        {
            for (int j = 0; j < width; j++)
            {
                x = da[i][j * 3];
                y = da[i][j * 3 + 1];
                z = da[i][j * 3 + 2];
                x = pow(x / 255, t) * 255;
                y = pow(y / 255, t) * 255;
                z = pow(z / 255, t) * 255;
                da[i][j * 3] = x;
                da[i][j * 3 + 1] = y;
                da[i][j * 3 + 2] = z;

            }
        }
}
void Matrix::nonmaxsupersion()
{
    double dp1 = 0, dp2 = 0, w = 0, x = 0, y = 0, g1 = 0, g2 = 0, g3 = 0, g4 = 0;
        for (int i = 1; i < height - 1; i++)
        {
            for (int j = 1; j < width - 1; j++)
            {
                x = (da[i][j * 3 + 3] - da[i + 1][3 * j] - da[i][3 * j] + da[i + 1][3 * j + 3]) / 2;
                y = (da[i][j * 3] - da[i + 1][3 * j] + da[i][3 * j + 3] - da[i + 1][3 * j + 3]) / 2;
                w = fabs(x) / fabs(y);
                if (fabs(x) < fabs(y))
                {
                    if (x * y > 0)
                    {
                        g1 = da[i - 1][(j - 1) * 3];
                        g2 = da[i - 1][j * 3];
                        g3 = da[i + 1][(j + 1) * 3];
                        g4 = da[i + 1][j * 3];
                    }
                    else
                    {
                        g1 = da[i - 1][(j + 1) * 3];
                        g2 = da[i - 1][j * 3];
                        g3 = da[i + 1][(j - 1) * 3];
                        g4 = da[i + 1][j * 3];
                    }
                }
                else
                {
                    if (x * y > 0)
                    {
                        g1 = da[i + 1][(j - 1) * 3];
                        g2 = da[i][(j - 1) * 3];
                        g3 = da[i - 1][(j + 1) * 3];
                        g4 = da[i][(j + 1) * 3];
                    }
                    else
                    {
                        g1 = da[i - 1][(j - 1) * 3];
                        g2 = da[i][(j - 1) * 3];
                        g3 = da[i + 1][(j + 1) * 3];
                        g4 = da[i][(j + 1) * 3];
                    }
                }
                dp1 = w * g1 + (1 - w) * g2;
                dp2 = w * g3 + (1 - w) * g4;
                if (da[i][j * 3] < dp1 || da[i][j * 3] < dp2)
                {
                    da[i][j * 3] = 0;
                    da[i][j * 3 + 1] = 0;
                    da[i][j * 3 + 2] = 0;
                }
            }
        }
}
void Matrix::doublefazhi(double lownumber, double heightnumber)
{
    if (lownumber > heightnumber)
        {
            return;
        }
        double x = 0;
        for (int i = 0; i < height; i++)
        {
            for (int j = 0; j < width; j++)
            {
                x = da[i][j * 3];
                if (x < lownumber)
                {
                    x = 0;
                }
                if (x <= heightnumber && x >= lownumber)
                {
                    for (int k = -1; k <= 3; k++)
                    {
                        for (int l = -1; l <= 3; l++)
                        {
                            if ((i + k) >= 0 && (i + k) < height && (j + l) >= 0 && (j + l) < width)
                            {
                                x = fmax(x, da[i + k][(j + l) * 3]);
                            }
                        }
                    }
                }
                da[i][j * 3] = x;
                da[i][j * 3 + 1] = x;
                da[i][j * 3 + 2] = x;
            }
        }
}
void Matrix::bianyuantiqu(double t, double lownumber, double heightnumber, int ksize, double strength)
{
    grey();
        gaosilvboqi(ksize, strength);
        tidu();
        nonmaxsupersion();
        if (t != 1)light(t);
        doublefazhi(lownumber, heightnumber);
}
void Matrix::ruihua(double t, double lownumber, double heightnumber, int ksize, double strength)
{
    Matrix c;
    c.copy(*this);
        c.bianyuantiqu(t, lownumber, heightnumber, ksize, strength);
        for (int i = 0; i < height; i++)
        {
            for (int j = 0; j < width * 3; j++)
            {
                da[i][j] += c.da[i][j];
                if (da[i][j] >= 256)
                {
                    da[i][j] = 255;
                }
            }
        }

}
void Matrix::junzhilvboqi(int ksize)
{
    if ((ksize % 2) == 0)return;
        int center = ksize / 2;
        int s = ksize * ksize;
        double sum = 0;
        double sum1 = 0, sum2 = 0;
        for (int i = 0; i < height; i++)
        {
            for (int j = 0; j < width; j++)
            {
                sum = 0;
                sum1 = 0;
                sum2 = 0;
                for (int k = i - center; k <= i + center; k++)
                {
                    for (int l = j - center; l <= j + center; l++)
                    {
                        if (l >= 0 && l < width && k >= 0 && k < height)
                        {
                            sum  += (da[k][l * 3] );
                            sum1 += (da[k][l * 3 + 1] );
                            sum2 += (da[k][l * 3 + 2] );
                        }
                    }
                }
                sum = sum / s;
                sum1 = sum1 / s;
                sum2 = sum2 / s;
                da[i][j * 3] = sum;
                da[i][j * 3 + 1] = sum1;
                da[i][j * 3 + 2] = sum2;
            }
        }
}
void Matrix::zhongzhilvboqi(int ksize)
{
    if (ksize % 2 == 0)return;
        int center = ksize / 2;
        double *sum=new double [ksize*ksize];
        double *sum1 = new double[ksize * ksize];
        double *sum2 = new double[ksize*ksize];
        int count = 0;
        for (int i = 0; i < height; i++)
        {
            for (int j = 0; j < width; j++)
            {
                for (int k = i - center; k <= i + center; k++)
                {
                    for (int l = j - center; l <= j + center; l++)
                    {
                        if (l >= 0 && l < width && k >= 0 && k < height)
                        {
                            sum [count]= (da[k][l * 3]);
                            sum1[count]= (da[k][l * 3 + 1]);
                            sum2[count]= (da[k][l * 3 + 2]);
                            count++;
                        }
                        else
                        {
                            sum[count] = 0;
                            sum1[count] = 0;
                            sum2[count] = 0;
                            count++;
                        }
                    }
                }
                count = 0;
                std::sort(sum, sum + ksize * ksize );
                std::sort(sum1, sum1 + ksize * ksize );
                std::sort(sum2, sum2 + ksize * ksize );
                da[i][j * 3] = sum[(ksize/2)*ksize+ksize/2];
                da[i][j * 3 + 1] = sum1[(ksize / 2) * ksize + ksize / 2];
                da[i][j * 3 + 2] = sum2[(ksize / 2) * ksize + ksize / 2];
            }
        }
}
